﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.LinearAlgebra;
using MentalAlchemy.Atomics;

namespace MentalAlchemy.Molecules
{
	/// <summary>
	/// [molecule]
	/// 
	/// Multivariate Gaussian distribution with arbitrary mean and covariance.
	/// </summary>
	public class GaussianDistribution
	{
		protected NormalDistribution stdDistr = new NormalDistribution(0, 1);
		protected float[,] sqrtEValues;
		protected float[,] vectorsValues;	// Q * L^{1/2}, Q -- matrix of eigenvectors, L -- matrix of eigenvalues.

		public bool FixCovariance = true;	// if provided covariance estimate turns out to be non-ositive definite then zero negative eigenvalues.

		public int Size { get; protected set; }
		public float[] Mean {get; protected set;}
		public float[,] Covariance{get; protected set;}
		public float[,] CovEigenvectors { get; protected set; }
		public float[,] CovEigenvalues { get; protected set; }


		#region - Public methods. -
		public GaussianDistribution(float[] mu, float[,] cov)
		{
			Size = mu.Length;
			Mean = mu;

			// compute eigenvectors and eigenvalues.
			SetCovariance(cov);
		}

		public float[] Next()
		{
			var res = new float[Size];
			for (int i = 0; i < Size; i++)
			{
				res[i] = (float)stdDistr.NextDouble();
			}
			res = MatrixMath.Mul(vectorsValues, res);
			res = VectorMath.Add(res, Mean);
			return res;
		} 

		public List<string> ToStrings ()
		{
			var res = new List<string>();
			res.Add("Mean:\t" + VectorMath.ConvertToString(Mean, '\t'));
			res.Add("Covariance matrix:");
			res.AddRange(MatrixMath.ConvertToRowsStringsList(Covariance, '\t'));
			res.Add("Eigenvalues:\t" + VectorMath.ConvertToString(MatrixMath.Diag(CovEigenvalues), '\t'));
			res.Add("Eigenvectors:");
			res.AddRange(MatrixMath.ConvertToRowsStringsList(MatrixMath.Transpose(CovEigenvectors), '\t'));
			return res;
		}
		#endregion

		#region - Protected methods. -
		protected void SetCovariance(float[,] cov)
		{
			Covariance = cov;
			var dcov = MatrixMath.ConvertToDoubles(cov);

			//var evDecomp = new EigenvalueDecomposition(Matrix.Create(dcov));
			//CovEigenvectors = MatrixMath.ConvertToFloats(evDecomp.EigenVectors.CopyToArray());
			//var evalues = evDecomp.EigenValues;
			//var realEvalues = VectorMath.Re(evalues);
			//CovEigenvalues = MatrixMath.Diagonal(realEvalues);

			//if (FixCovariance)
			//{
			//    for (int i = 0; i < cov.GetLength(0); i++)
			//    {
			//        if (CovEigenvalues[i, i] < 0) CovEigenvalues[i, i] = 0;
			//    }
			//}

			//sqrtEValues = MatrixMath.Sqrt(CovEigenvalues);
			//vectorsValues = MatrixMath.Mul(CovEigenvectors, sqrtEValues);

			var eigv = new double[0];
			var eigvec = new double[0,0];
			alglib.evd.smatrixevd(dcov, dcov.GetLength(0), 1, false, ref eigv, ref eigvec);

			var evals = VectorMath.ConvertToFloats(eigv);
			CovEigenvalues = MatrixMath.Diagonal(evals);
			CovEigenvectors = MatrixMath.ConvertToFloats(eigvec);

			if (FixCovariance)
			{
				for (int i = 0; i < cov.GetLength(0); i++)
				{
					if (CovEigenvalues[i, i] < 0) CovEigenvalues[i, i] = -CovEigenvalues[i, i];
				}
			}

			sqrtEValues = MatrixMath.Sqrt(CovEigenvalues);
			vectorsValues = MatrixMath.Mul(CovEigenvectors, sqrtEValues);

		} 
		#endregion
	}
}
